home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / MARSAG.C < prev    next >
C/C++ Source or Header  |  1996-03-26  |  9KB  |  284 lines
  1. /* marsag.c */
  2. /* -------------------------------------------------------------------- *
  3.  * This random number generator originally appeared in "Toward a        *
  4.  *   Universal Random Number Generator" by George Marsaglia and Arif    *
  5.  *   Zaman. Florida State University Report: FSU-SCRI-87-50 (1987)    *
  6.  *                                    *
  7.  * It was later modified by F. James and published in "A Review of      *
  8.  *   Pseudo-random Number Generators"                                   *
  9.  *                                    *
  10.  * THIS IS THE BEST KNOWN RANDOM NUMBER GENERATOR AVAILABLE.        *
  11.  *   (But, a newly discovered technique can yield a period of 10^600.    *
  12.  *    But that is still in the development stage.)            *
  13.  *                                    *
  14.  * It passes ALL of the tests for random number generators and has a    *
  15.  *   period of 2^144, is portable (gives bit-identical results on all    *
  16.  *   machines with at least 24-bit mantissas in the floating point    *
  17.  *   representation).                            *
  18.  *                                    *
  19.  * The algorithm is a combination of a Fibonacci sequence (with lags of *
  20.  *   97 and 33, an operation "subtraction plus one, modulo one") and an *
  21.  *   "arithmetic sequence" (using subtraction).             *
  22.  * ==================================================================== *
  23.  * This C language version was written by Jim Butler, and is based on a *
  24.  * FORTRAN program posted by David LaSalle of Florida State University. *
  25.  * -------------------------------------------------------------------- */
  26. #include <stdlib.h>
  27.  
  28. #define FULL_SIZE ((unsigned)RAND_MAX + 1U)
  29.  
  30. /* ------------------- */
  31. /* FUNCTION PROTOTYPES */
  32. /* ------------------- */
  33. # undef F
  34. # if defined(__STDC__) || defined(__PROTO__)
  35. #    define    F( P )    P
  36. # else
  37. #    define    F( P )    ()
  38. # endif
  39. /* INDENT OFF */
  40. extern    double    Dranmar F((void));
  41. extern    int    Iranmar F((void));
  42. extern    void    ranmar F((double *, int));
  43. extern    void    rmarin F((unsigned, unsigned));
  44.  
  45. # undef F
  46. /* INDENT ON */
  47. # if defined(TEST_MARSAG)
  48. #include <stdio.h>
  49.  
  50. void
  51. main()
  52. {
  53.     double  temp[100];
  54.     int     ij, kl, len, i;
  55.  
  56.     /* random seeds for the test case: */
  57.     ij = 1802;
  58.     kl = 9373;
  59. #if 0
  60.     /* do the initialization */
  61.     rmarin(ij, kl);
  62. # endif
  63.     /* generate 20,000 random numbers */
  64.     len = 100;
  65.     for (i = 1; i <= 200; i++)
  66.     ranmar(temp, len);
  67.  
  68. /*
  69.  * If the random number generator is working properly, the next six
  70.  * random numbers should be:
  71.  *
  72.  *  6533892.0  14220222.0  7275067.0  6172232.0  8354498.0  10633180.0
  73.  */
  74.     len = 6;
  75.     ranmar(temp, len);
  76.  
  77.     for (i = 0; i < 6; i++)
  78.     printf("%12.1f", 4096.0 * 4096.0 * temp[i]);
  79.  
  80.     printf("\nCompare:\n   6533892.0  14220222.0   7275067.0 "
  81.     "  6172232.0   8354498.0  10633180.0\n");
  82. }
  83. # endif
  84. /* ---------------------------------------------------------- */
  85. /* Initialized Table.  These numbers have values identical to */
  86. /* those created when function rmarin() is called as follows: */
  87. /*                                  */
  88. /*-    int    ij, kl;                       */
  89. /*                                  */
  90. /*-    ij = 1802;                          */
  91. /*-    kl = 9373;                          */
  92. /*                                  */
  93. /*    rmarin(ij, kl);                       */
  94. /* ---------------------------------------------------------- */
  95. /* INDENT OFF */
  96.  
  97. #define _D_    16777216.0
  98. static    double    u[98] =
  99. {
  100.      0.0      ,  13697435.0 / _D_,       3833429.0 / _D_,
  101.   12353926.0 / _D_,   2287754.0 / _D_,       3468638.0 / _D_,
  102.    1232959.0 / _D_,   8059805.0 / _D_,      10745739.0 / _D_,
  103.    4236676.0 / _D_,   2095136.0 / _D_,       1349346.0 / _D_,
  104.    3672867.0 / _D_,  14563641.0 / _D_,      15473517.0 / _D_,
  105.    9897259.0 / _D_,   2207061.0 / _D_,        929657.0 / _D_,
  106.    8109095.0 / _D_,   5246947.0 / _D_,       1066111.0 / _D_,
  107.    8460236.0 / _D_,  13162386.0 / _D_,        501474.0 / _D_,
  108.   10402355.0 / _D_,    352505.0 / _D_,       2104170.0 / _D_,
  109.   12045925.0 / _D_,   4350943.0 / _D_,      13996856.0 / _D_,
  110.    9897761.0 / _D_,   6626452.0 / _D_,      15057436.0 / _D_,
  111.    3168599.0 / _D_,  14038489.0 / _D_,       8550848.0 / _D_,
  112.    5242835.0 / _D_,  13296102.0 / _D_,      11969002.0 / _D_,
  113.      95246.0 / _D_,   5917978.0 / _D_,       8555838.0 / _D_,
  114.   13557738.0 / _D_,   1526088.0 / _D_,      11197237.0 / _D_,
  115.   15721125.0 / _D_,  14247931.0 / _D_,        897046.0 / _D_,
  116.   15537441.0 / _D_,  16645456.0 / _D_,      16279884.0 / _D_,
  117.    1289925.0 / _D_,  14032128.0 / _D_,      10641039.0 / _D_,
  118.    9961793.0 / _D_,   2737638.0 / _D_,       5073398.0 / _D_,
  119.    5231619.0 / _D_,   2007688.0 / _D_,      15753584.0 / _D_,
  120.   12368695.0 / _D_,  12926325.0 / _D_,      10522018.0 / _D_,
  121.    8692194.0 / _D_,   8531802.0 / _D_,      14755384.0 / _D_,
  122.     276334.0 / _D_,   9157821.0 / _D_,        989353.0 / _D_,
  123.    6093627.0 / _D_,  15866666.0 / _D_,       9532882.0 / _D_,
  124.    3434034.0 / _D_,    710155.0 / _D_,        672726.0 / _D_,
  125.   12734991.0 / _D_,  13809842.0 / _D_,       4832132.0 / _D_,
  126.    9753458.0 / _D_,  11325486.0 / _D_,      12137466.0 / _D_,
  127.    3617374.0 / _D_,   4913050.0 / _D_,       9978642.0 / _D_,
  128.   12740205.0 / _D_,  15754026.0 / _D_,       4928136.0 / _D_,
  129.    8545553.0 / _D_,  12893795.0 / _D_,       8164497.0 / _D_,
  130.   12420478.0 / _D_,   8192378.0 / _D_,       2028808.0 / _D_,
  131.    1183983.0 / _D_,   3474722.0 / _D_,      15616920.0 / _D_,
  132.   16298670.0 / _D_,  14606645.0 / _D_
  133. };
  134.  
  135. static    double    c  =   362436.0 / 16777216.0,
  136.         cd =  7654321.0 / 16777216.0,
  137.         cm = 16777213.0 / 16777216.0;
  138.  
  139. static    int    i97 = 97, j97 = 33;
  140.  
  141. /* INDENT ON */
  142.  
  143. /* -------------------------------------------------------------------- *
  144.  * This is the initialization routine for the random number generator    *
  145.  * RANMAR(). NOTE: The seed variables can have values between:        *
  146.  *                                    *
  147.  *    0 <= IJ <= 31328                            *
  148.  *    0 <= KL <= 30081                            *
  149.  *                                    *
  150.  * The random number sequences created by these two seeds are long    *
  151.  * enough to complete an entire calculation with. For example, if    *
  152.  * several different groups are working on different parts of the    *
  153.  * same calculation, each group could be assigned its own IJ seed.    *
  154.  * This would leave each group with 30000 choices for the second    *
  155.  * seed. That is to say, this random number generator can create    *
  156.  * 900 million different subsequences, each subsequence having a    *
  157.  * length of approximately 10^30.                    *
  158.  *                                    *
  159.  * Use IJ = 1802 & KL = 9373 to test the random number generator.    *
  160.  * First, subroutine RANMAR should be used to generate 20000 random    *
  161.  * numbers. Then display the next six random numbers generated each    *
  162.  * multiplied by 4096*4096. If the random number generator is working    *
  163.  * properly, the random numbers should be:                *
  164.  *                                    *
  165.  *  6533892.0  14220222.0  7275067.0  6172232.0  8354498.0  10633180.0    *
  166.  * -------------------------------------------------------------------- */
  167. # if defined(__STDC__) || defined(__PROTO__)
  168. void
  169. rmarin(unsigned IJ, unsigned KL)
  170. # else
  171. void
  172. rmarin(IJ, KL)
  173. unsigned IJ;
  174. unsigned KL;
  175. # endif
  176. {
  177.     int     i, ij, j, k, kl, l, ii, jj, m;
  178.     double  s, t;
  179.  
  180.     ij = IJ % 31329;
  181.     kl = KL % 30082;
  182.  
  183.     i = (ij / 177) % 177 + 2;
  184.     j = ij % 177 + 2;
  185.     k = (kl / 169) % 178 + 1;
  186.     l = kl % 169;
  187.  
  188.     for (ii = 1; ii <= 97; ii++)
  189.     {
  190.     s = 0.0;
  191.     t = 0.5;
  192.     for (jj = 1; jj <= 24; jj++)
  193.     {
  194.         m = (((i * j) % 179) * k) % 179;
  195.         i = j;
  196.         j = k;
  197.         k = m;
  198.         l = (53 * l + 1) % 169;
  199.         if ((l * m) % 64 >= 32)
  200.         s += t;
  201.         t *= 0.5;
  202.     }
  203.     u[ii] = s;
  204.     }
  205.                     /* INDENT OFF */
  206.     c  =   362436.0 / 16777216.0;
  207.     cd =  7654321.0 / 16777216.0;
  208.     cm = 16777213.0 / 16777216.0;
  209.  
  210.                     /* INDENT ON */
  211.     i97 = 97;
  212.     j97 = 33;
  213. }
  214. /* ==================================================================== */
  215. /* Dranmar - returns pseudo random double in [0, 1)            */
  216. /* ==================================================================== */
  217. double
  218. Dranmar()
  219. {
  220.     double  uni;
  221.  
  222.     uni = u[i97] - u[j97];
  223.     if (uni < 0.0)
  224.     {
  225.     uni += 1.0;
  226.     }
  227.     u[i97] = uni;
  228.     i97--;
  229.     if (i97 == 0)
  230.     {
  231.     i97 = 97;
  232.     }
  233.     j97--;
  234.     if (j97 == 0)
  235.     {
  236.     j97 = 97;
  237.     }
  238.     c -= cd;
  239.     if (c < 0.0)
  240.     {
  241.     c += cm;
  242.     }
  243.     uni -= c;
  244.     if (uni < 0.0)
  245.     {
  246.     uni += 1.0;
  247.     }
  248.     return uni;
  249. }
  250. /* ==================================================================== */
  251. /* Iranmar - returns pseudo random integer in [0, RAND_MAX]        */
  252. /* ==================================================================== */
  253. int
  254. Iranmar()
  255. {
  256.     return((int) ((double)FULL_SIZE * Dranmar()) & RAND_MAX);
  257. }
  258. /* -------------------------------------------- */
  259. /* ranmar - places len random numbers in rvec[] */
  260. /* -------------------------------------------- */
  261. # if defined(__STDC__) || defined(__PROTO__)
  262. void
  263. ranmar(double *rvec, int len)
  264. # else
  265. void
  266. ranmar(rvec, len)
  267. double *rvec;
  268. int    len;
  269. # endif
  270. /* -------------------------------------------------------------------- *
  271.  * This is the random number generator proposed by George Marsaglia    *
  272.  * in Florida State University Report: FSU-SCRI-87-50.            *
  273.  * It was slightly modified by F. James to produce an array of        *
  274.  * pseudorandom numbers.                        *
  275.  * -------------------------------------------------------------------- */
  276. {
  277.     int     ivec;
  278.  
  279.     for (ivec = 0; ivec < len; ++ivec)
  280.     {
  281.     rvec[ivec] = Dranmar();
  282.     }
  283. }
  284.